Modello Statistico per la Previsione del Peso Neonatale

1 Introduzione

L’analisi nasce dall’esigenza di avere una previsione anticipata del peso del neonato, al fine di una più corretta pianificazione delle risorse da allocare in caso di peso inferiore alle aspettative.

Il peso è difatti una caratteristica indicativa e di grande rilievo per la salute del neonato e in caso sia troppo basso, devono essere predisposte strumentazioni e persone per effettuare attività ben precise, atte a riportare il peso entro i parametri corretti.

Il progetto si inserisce all’interno di un contesto di crescente attenzione verso la prevenzione delle complicazioni neonatali. La possibilità di prevedere il peso alla nascita dei neonati rappresenta un’opportunità fondamentale per migliorare la pianificazione clinica e ridurre i rischi associati a nascite problematiche, come parti prematuri o neonati con basso peso. Di seguito, i principali benefici che questo progetto porterà all’azienda e al settore sanitario:

  1. Miglioramento delle previsioni cliniche:

    • Il peso del neonato è un indicatore chiave della sua salute. Avere un modello predittivo accurato consente al personale medico di intervenire tempestivamente in caso di anomalie, riducendo le complicazioni perinatali come le difficoltà respiratorie o l’ipoglicemia.
  2. Ottimizzazione delle risorse ospedaliere:

    • Sapere in anticipo quali neonati potrebbero avere bisogno di cure intensive aiuta a organizzare le risorse umane e tecnologiche degli ospedali in modo efficiente. Questo si traduce in una riduzione dei costi operativi e una migliore pianificazione dell’utilizzo delle unità di terapia intensiva neonatale (TIN).
  3. Prevenzione e identificazione dei fattori di rischio:

    • Il modello potrà evidenziare i fattori che maggiormente influenzano negativamente il peso del neonato (come il fumo materno, gravidanze multiple o età avanzata della madre). Queste informazioni sono preziose per la prevenzione e la gestione personalizzata delle gravidanze, permettendo interventi proattivi in caso di rischio elevato.
  4. Valutazione delle pratiche ospedaliere:

    • Attraverso un’analisi comparativa tra i tre ospedali coinvolti, l’azienda potrà identificare eventuali differenze nei risultati clinici, come una maggiore incidenza di parti cesarei in una determinata struttura. Ciò consente di monitorare la qualità delle pratiche e armonizzare i protocolli tra i diversi centri ospedalieri, migliorando la coerenza delle cure.
  5. Supporto alla pianificazione strategica:

    • L’analisi dei dati e le previsioni possono essere utilizzate per prendere decisioni informate non solo a livello clinico ma anche strategico. L’azienda potrà sfruttare queste informazioni per implementare nuove politiche di salute pubblica, garantendo un impatto positivo sui tassi di mortalità e morbilità neonatale.

2 Analisi descrittiva

In questa sezione si utilizzeranno i più noti indicatori statistici per descrivere il dataset utilizzato al fine di:

  • conoscerne adeguatamente le variabili,

  • individuare eventuali valori anomali e inquadrarli

  • comprendere la struttura generale dei dati

2.1 Il dataset

setwd("~/ProfessionAI - Data Science/03_Statistica inferenziale")
neonati = read.csv("neonati.csv", stringsAsFactors = T)

2.2 Le librerie

Di seguito le librerie utilizzate per condurre l’analisi. Sostanzialmente tutto il pacchetto tidyverse per analisi, manipolazione dei dati, nonchè creazioni di grafici.

Le librerie car, MASS e lmtest tutte dedicate ai test sul modello statistico definito.

library(ggplot2)
library(dplyr)
library(patchwork)
library(moments)
library(tidyr)
library(car)
library(MASS)
library(lmtest)

2.3 Le variabili

I dati analizzati sono stati raccolti da 3 ospedali differenti, in ciascuno dei quali le variabili registrate sono state 8, con 2500 osservazioni totali.

Queste sono le 8 variabili studiate all’interno del campione

colnames(neonati)
 [1] "Anni.madre"   "N.gravidanze" "Fumatrici"    "Gestazione"   "Peso"        
 [6] "Lunghezza"    "Cranio"       "Tipo.parto"   "Ospedale"     "Sesso"       

Di seguito una loro breve descrizione:

  • Anni.madre: Misura dell’età in anni.

  • N.gravidanze: Quante gravidanze ha avuto la madre.

  • Fumatrici: Un indicatore binario (0=non fumatrice, 1=fumatrice).

  • Gestazione: Numero di settimane di gestazione.

  • Peso: Peso alla nascita in grammi.

  • Lunghezza: lunghezza del neonato in mm, misurabile anche durante la gravidanza tramite ecografie.

  • Cranio: diametro craniale, misurabile anche durante la gravidanza tramite ecografie.

  • Tipo.parto: Naturale o cesareo.

  • Ospedale: Ospedale di nascita 1, 2 o 3.

  • Sesso: Maschio (M) o femmina (F).

Di seguito invece una visione di sintesi dei quantili e della media di ciascuna variabile nel dataset

summary(neonati)
   Anni.madre     N.gravidanze       Fumatrici        Gestazione   
 Min.   : 0.00   Min.   : 0.0000   Min.   :0.0000   Min.   :25.00  
 1st Qu.:25.00   1st Qu.: 0.0000   1st Qu.:0.0000   1st Qu.:38.00  
 Median :28.00   Median : 1.0000   Median :0.0000   Median :39.00  
 Mean   :28.16   Mean   : 0.9812   Mean   :0.0416   Mean   :38.98  
 3rd Qu.:32.00   3rd Qu.: 1.0000   3rd Qu.:0.0000   3rd Qu.:40.00  
 Max.   :46.00   Max.   :12.0000   Max.   :1.0000   Max.   :43.00  
      Peso        Lunghezza         Cranio    Tipo.parto Ospedale   Sesso   
 Min.   : 830   Min.   :310.0   Min.   :235   Ces: 728   osp1:816   F:1256  
 1st Qu.:2990   1st Qu.:480.0   1st Qu.:330   Nat:1772   osp2:849   M:1244  
 Median :3300   Median :500.0   Median :340              osp3:835           
 Mean   :3284   Mean   :494.7   Mean   :340                                 
 3rd Qu.:3620   3rd Qu.:510.0   3rd Qu.:350                                 
 Max.   :4930   Max.   :565.0   Max.   :390                                 

Balza subito all’occhio che ci sia qualcosa che non vada nella variabile Anni.madre, poichè presenta un valore minimo di 0, cosa ovviamente impossibile.

Si effettua quindi una verifica più puntuale per individuare altri valori chiaramente anonimi, filtrando il dataset per la variabile anni.madre < 15

neonati %>%
  filter(Anni.madre < 15) %>%
  arrange(Anni.madre)
Anni.madre N.gravidanze Fumatrici Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
0 0 0 39 3060 490 330 Nat osp3 M
1 1 0 41 3250 490 350 Nat osp2 F
13 0 0 38 2760 470 325 Nat osp2 F
14 1 0 39 3510 490 365 Nat osp2 M
14 0 0 39 3550 500 355 Ces osp1 M

Si evince che, pur essendo valori di età bassi, 13 e 14 non possono essere scartati perchè comunque verosimili. Lo stesso non si può ovviamente dire per 1 e 0.

Probabilmente si tratta di errori di battitura che però potrebbero influenzare negativamente le conclusioni che si prenderanno. Si decide dunque di eliminare queste 2 osservazioni dal dataset

neonati.clean = neonati %>%
  filter(Anni.madre > 1)

2.4 Gli indici

Di seguito verranno calcolati gli indici di variabilità e di forma, così da avere un’idea della possibile distribuzione delle variabili e soffermarci su envtuali valori degni di nota.

Resteranno fuori da questa analisi quantitativa, le variabili qualitative come Sesso e Ospedale.

2.4.1 Indici di variabilità

Gli indici di variabilità forniscono un secondo livello di informazioni sulle variabili, in particolare indicano quanto esse siano disperse o concentrate attorno al valore medio.

Di seguito una tabella che raccoglie, per ogni variabile quantitativa, gli indici di variabilità corrispettivi.

attach(neonati.clean)

variance = round(c(var(Anni.madre),
var(N.gravidanze),
var(Gestazione),
var(Peso),
var(Lunghezza),
var(Cranio)),2)

standard_dev = round(sqrt(variance),2)

IQR = round(c(IQR(Anni.madre),
IQR(N.gravidanze),
IQR(Gestazione),
IQR(Peso),
IQR(Lunghezza),
IQR(Cranio)),2)

variation_index = round(c(sd(Anni.madre)/mean(Anni.madre),
                    sd(N.gravidanze)/mean(N.gravidanze),
                    sd(Gestazione)/mean(Gestazione),
                    sd(Peso)/mean(Peso),
                    sd(Lunghezza)/mean(Lunghezza),
                    sd(Cranio)/mean(Cranio)),2)

VARIATION_GLOBAL = data.frame(IQR, variance, standard_dev, variation_index)
rows = c("anni_madre", "n_gravidanze", "gestazione", "peso", "lunghezza", "cranio")
row.names(VARIATION_GLOBAL) = rows

VARIATION_GLOBAL
Table 1: Indici di variazione
IQR variance standard_dev variation_index
anni_madre 7 27.22 5.22 0.19
n_gravidanze 1 1.64 1.28 1.30
gestazione 2 3.49 1.87 0.05
peso 630 275865.90 525.23 0.16
lunghezza 30 693.21 26.33 0.05
cranio 20 269.93 16.43 0.05

2.4.2 Indici di forma

Gli indici di forma forniscono informazioni sulla forma della distribuzione, indicando quanto essa sia asimmetrica, con il vertice della campana dunque spostato verso sinistra o destra rispetto a una gaussiana standard, oppure quanto tale vertice si alzi o abbassi rispetto a una gaussiana standard.

Di seguito una tabella con l’indice di asimmetria e curtosi per ogni variabile

Skewness = round(c(skewness(Anni.madre),
  skewness(N.gravidanze),
  skewness(Gestazione),
  skewness(Peso),
  skewness(Lunghezza),
  skewness(Cranio)),2)

Kurtosis = round(c(kurtosis(Anni.madre),
             kurtosis(N.gravidanze),
             kurtosis(Gestazione),
             kurtosis(Peso),
             kurtosis(Lunghezza),
             kurtosis(Cranio)),2)

SHAPE_GLOBAL = data.frame(Skewness,Kurtosis)
row.names(SHAPE_GLOBAL) = rows

SHAPE_GLOBAL
Table 2: Indici di forma
Skewness Kurtosis
anni_madre 0.15 2.89
n_gravidanze 2.51 13.98
gestazione -2.07 11.26
peso -0.65 5.03
lunghezza -1.51 9.48
cranio -0.79 5.94

2.4.3 Commenti sugli indici

Possiamo senza dubbio concludere che le variabili osservate non presentano un’eccessiva varibilità. Confrontando l’indice di variazione infatti, vale a dire la media divisa per la deviazione standard, solo il numero di gravidanze sembra mostrare un particolare scostamento, con la deviazione standard il 30% più grande della media.

Il numero di gravidanze, insieme alle settimane di gestazione, presentano invece la maggiore asimmetria. In particolare:

  • il numero di gravidanze ha un’asimmetria positiva, verso sinistra dunque. Questo ci dice che, tendenzialmente, il nostro campione è composto da mamme con pochi parti alle spalle. La coda destra è comunque lunga, ad indicare che, anche se in poche occorrenze, c’è chi era alla 13ma gestazione.
Figure 1
  • Il numero di settimane di gestazione invece ha asimmetria negativa, vale a dire verso destra. Il risultato non stupisce, poichè il numero tipico di settimane di gestazione è 40. Essendo un valore “massimo”, non può generare una gaussiana standard, altrimenti si avrebbe il 50% di probabilità di avere parti oltre la 40esima. Va posta attenzione sul fatto che però ci sia una coda sinistra pronunciata, con casi di nascita addirittura alla 25ma settimana. Questo potrebbe complicare le cose in fase di modellazione.
Figure 2

La Kurtosi invece ci indica la chiara tendenza di tutte le variabili a essere lepticurtiche, cioè con una “campana” più stretta e alta rispetto a una normale standard. Proprio come per l’asimmetria, le variabili decisamente più lepticurtiche, sono proprio le settimane di gestazione e il numero di gravidanze pregresse.

3 Test statistici sul campione

In questo paragrafo verranno sondate alcune ipotesi che i dati lasciano supporre. Questo passo si rende necessario perchè, nonostante i dati del campione possano puntare in una certa direzione, bisogna sempre ricordare che, anche se numeroso, si tratta sempre di un campione.

Pertanto, prima di generalizzare quello che il campione mostra, si deve cercare di capire se ci sia margine o meno per effettuare inferenza.

3.1 Differenze fisiche tra i sessi

La prima domanda cui vogliamo dare risposta è: esiste differenza nelle misure antropometriche (lunghezza, peso e diametro craniale) tra i due sessi?

Si pone nuovamente l’attenzione sul fatto che la domanda sia generale. Non si sta chiedendo se esistano nel campione, delle differenze tra sessi. Per quello un grafico sarebbe più che sufficiente.

Quello che si sta chiedendo è se si possa estendere a una popolazione quanto visto sul campione.

Innanzitutto, vediamo quale sia l’andamento nel campione, relativamente alle misure antropometriche, per i due sessi.

Figure 3
Figure 4
Figure 5

I 3 grafici mostrano che, nel campione, esiste una differenza tra maschi e femmine, per tutte le misure antropometriche, laddove i maschi le hanno tendenzialmente maggiori rispetto alle femmine.

Degno di nota il fatto che esistono molti valori non racchiusi dal boxplot e che dovranno essere studiati in fase di definizione del modello lineare.

Testiamo ora il seguente sistema di ipotesi:

\[ \begin{cases} H_{0} : \mu_{M} - \mu_{F} = 0, & \mbox{La media per i maschi } \mu_{M}\mbox{ coincide coincide con la media per le femmine } \mu_{F}\\ H_{1} : \mu_{M} - \mu_{F} \neq 0, & \mbox{Le medie tra maschi e femmine sono diverse } \end{cases} \]

Di volta in volta assumeremo la media per i 3 indicatori antropometrici: lunghezza, peso e diametro craniale.

Si userà il test t, essendo il più adatto per confronti di medie tra due gruppi. Il livello di significatività (cioè la probabilità di rifiutare l’ipotesi nulla quando in realtà è vera) sarà fissato al 5%.

3.1.1 Differenza di peso tra i sessi

Il test verificherà che la media di peso dei maschi, sottratta alla media di peso delle femmine sia uguale a 0, vale a dire che le due media siano identiche.

Il test assume che le due medie siano uguali. Se questo fosse vero, la differenza tra le due, estraendo dei campioni, si distribuirebbe secondo una curva t di Student.

t.test(Peso~Sesso, data = neonati.clean) 

    Welch Two Sample t-test

data:  Peso by Sesso
t = -12.115, df = 2488.7, p-value < 2.2e-16
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
 -287.4841 -207.3844
sample estimates:
mean in group F mean in group M 
       3161.061        3408.496 

Il test mostra un p-value praticamente a 0. Questo valore indica la probabilità di ottenere una differenza tra le medie uguale rispetto a quella vista in questo campione o peggio, se fosse vera l’ipotesi nulla.

Questo vuol dire che la probabilità che il nostro campione presenti una differenza tra le medie così diversa da 0 rasenta l’impossibilità. Talmente improbabile che diventa assolutamente più probabile che sia sbagliata l’ipotesi per cui le 2 medie sono coincidenti.

Ergo, il test dimostra che le due medie sono diverse tra loro in maniera statisticamente molto significativa.

Possiamo pertanto estendere a tutta la popolazione l’affermazione che: c’è differenza di peso tra i neonati a dipendere dal sesso.

3.1.2 Differenza di diametro craniale tra i sessi

t.test(Cranio~Sesso, data = neonati.clean)  

    Welch Two Sample t-test

data:  Cranio by Sesso
t = -7.4366, df = 2489.4, p-value = 1.414e-13
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
 -6.110504 -3.560417
sample estimates:
mean in group F mean in group M 
       337.6231        342.4586 

Anche in questo caso abbiamo un p-value vicinissimo allo 0. Si rifiuta pertanto l’ipotesi nulla di uguaglianza tra le due medie e possiamo concludere che esiste differenza, tra maschi e femmine, nel diametro craniale.

3.1.3 Differenza di lunghezza tra i sessi

t.test(Lunghezza~Sesso, data = neonati.clean)              

    Welch Two Sample t-test

data:  Lunghezza by Sesso
t = -9.5823, df = 2457.3, p-value < 2.2e-16
alternative hypothesis: true difference in means between group F and group M is not equal to 0
95 percent confidence interval:
 -11.939001  -7.882672
sample estimates:
mean in group F mean in group M 
       489.7641        499.6750 

Anche in questo caso il p-value è vicino allo 0 e si rifiuta l’ipotesi nulla. Anche per la lunghezza possiamo affermare che esiste differenza tra maschi e femmine in generale.

3.2 Peso e lunghezza del campione vs popolazione

Vogliamo verificare che il dato medio di peso e lunghezza del nostro campione, sia allineato ai valori medi della popolazione. Il livello di significatività lo si pone pari al 5%.

\[ \begin{cases} H_{0} : \mbox{La media del campione coincide con la media della popolazione } \\ H_{1} : \mbox{La media del campione è diversa dalla media della popolazione } \end{cases} \]

Per farlo si utilizzerà sempre un test t, perchè non conoscendo la deviazione standard della popolazione, la curva t di Student userà quella del campione e incorporerà l’incertezza di questa approssimazione.

Per ottenere il dato medio di peso e lunghezza dei neonati, che sarà il nostro dato di popolazione, si è attinto al sito dell’ospedale pediatrico Bambino Gesù per cui risulta:

  • peso medio: 3300 grammi

  • lunghezza media: 50 cm

t.test(Peso, alternative = "two.sided", mu = 3300)

    One Sample t-test

data:  Peso
t = -1.505, df = 2497, p-value = 0.1324
alternative hypothesis: true mean is not equal to 3300
95 percent confidence interval:
 3263.577 3304.791
sample estimates:
mean of x 
 3284.184 

Il p-value supera di poco il 13% ed essendo quindi maggiore del livello di significatività, possiamo accettare l’ipotesi nulla: la media del campione è allineata con quella della popolazione

t.test(Lunghezza, alternative = "two.sided", mu = 500)

    One Sample t-test

data:  Lunghezza
t = -10.069, df = 2497, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 500
95 percent confidence interval:
 493.6628 495.7287
sample estimates:
mean of x 
 494.6958 

Lo stesso non può essere detto della lunghezza, per cui il p-value è molto vicino allo 0 e dobbiamo pertanto rifiutare l’ipotesi nulla: le lunghezze di questo campione non sono allineate con la media della popolazione.

3.3 Parti naturali vs cesarei

Concentriamoci innanzitutto sul descrivere queste variabili all’interno del nostro campione.

Tipo.parto
 Ces  Nat 
 728 1770 

Il dataset mostra una netta maggioranza di parti naturali che diventa ancora più ovvia con l’ausilio di un grafico a barre

Figure 6

Provando ad ampliare la vista includendo ora anche gli ospedali abbiamo questo grafico:

Figure 7

Da questo grafico capiamo che nel campione sembra che gli ospedali seguano lo stesso andamento sui parti, quasi fossero equidistribuiti sui 3. Verrà saggiata questa ipotesi con dei test statistici ad hoc.

3.3.1 Si fanno più parti naturali che cesarei?

Il grafico sembrerebbe assolutamente indicare di si, ma possiamo essere così certi solo per quanto riguarda il nostro campione. Se volessimo generalizzare questa affermazione dovremmo verificare le probabilità di aver pescato un campione come il nostro.

Il test più adeguato per i confronti tra variabili categoriche è il Chi-quadrato. Può essere usato per verificare l’indipendenza di due variabili, o per verificare se una variabile si scosti significativamente da una distribuzione attesa (goodness of fit).

È proprio con questa seconda accezione che effettueremo il test, assumendo che non si facciano più parti di un tipo che di un altro. Quindi:

\[ \begin{cases} H_{0} : P\left( Cesareo \right) = 0,5 & \\ H_{1} : P\left( Cesareo \right) \neq 0,5\end{cases} \]

chisq.test(table(Tipo.parto))

    Chi-squared test for given probabilities

data:  table(Tipo.parto)
X-squared = 434.65, df = 1, p-value < 2.2e-16

Il p-value è molto vicino allo 0, dobbiamo quindi rifiutare l’ipotesi nulla e possiamo affermare che effettivamente sussiste uno sbilanciamento nei tipi di parto, a favore dei naturali piuttosto che dei cesarei.

3.3.2 In alcuni ospedali si fanno più cesarei?

Dal momento che abbiamo dimostrato che non si fanno più cesarei che naturali, la domanda è da intendersi: ” C’è un ospedale che fa più cesarei degli altri? ”

Il grafico Figure 7 mostra che tra gli ospedali sembra esserci un’equidistribuzione dei cesarei, ma lasceremo al test saggiare questa ipotesi.

PartixOsp = neonati.clean %>%
  dplyr::select(Ospedale, Tipo.parto)

PartixOsp %>%
  filter(Tipo.parto == "Ces") %>%
  mutate(Tipo.parto = droplevels(Tipo.parto)) %>%
  table() %>%
  chisq.test()

    Chi-squared test for given probabilities

data:  .
X-squared = 1, df = 2, p-value = 0.6065

Il p-value è molto alto, poco oltre il 60%. Non rifiutiamo pertanto l’ipotesi nulla e affermiamo che tra gli ospedali c’è una distribuzione equa dei cesarei. Nessun ospedale ne fa quindi più degli altri in maniera statisticamente significativa.

3.3.3 Indipendenza delle variabili Ospedale e Tipo.parto

Infine, verifichiamo se le due variabili Ospedale e Tipo.parto siano indipendenti. La loro indipendenza implicherebbe che, conoscere il valore di una non consentirebbe di ricavare il valore dell’altra.

Il sistema di ipotesi pertanto è:

\[ \begin{cases} H_{0} : \mbox {Tipo.parto e Ospedale sono indipendenti} & \\ H_{1} : \mbox{Tipo.parto e Ospedale non sono indipendenti}\end{cases} \]

PartixOsp %>%
  table() %>%
  chisq.test()

    Pearson's Chi-squared test

data:  .
X-squared = 1.083, df = 2, p-value = 0.5819

Il p-value è quasi al 60%, molto alto. Possiamo quindi accettare l’ipotesi nulla e affermare che le variabili Tipo.parto e Ospedale sono indipendenti. Questo vuol dire che non osserviamo nessuno schema ricorrente per cui, conosciuto l’ospedale, possiamo dire quale tipo di parto sarà praticato e viceversa.

4 Il modello di regressione

In questo capitolo verrà illustrato il procedimento seguito per definire un modello lineare basato sulle variabili del dataset e che consenta di prevedere in anticipo il peso del neonato.

Verrà seguita la procedura iterativa stepwise, partendo dal modello più complesso, con un regressore per ogni variabile del dataset e da lì si procederà per scremature sequenziali fino ad ottenere un buon compromesso tra varianza e semplicità.

I test e i criteri di valutazione che accompagneranno la procedura saranno:

  • il t test sui coefficienti di regressione delle variabili e relativo p-value. Il sistema di ipotesi è il seguente

\[ \begin{cases} H_{0} : \beta_{i} = 0 &\\ H_{0} : \beta_{i} \neq 0 \end{cases} \]

  • Il coefficiente di determinazione aggiustato R^2 di ogni modello che indica la percentuale di adattabilità del modello ai dati veri

  • Il criterio di informazione Bayesiano (BIC)

  • Il test ANOVA

4.1 Verifiche preliminari

Il dataset si presenta con la variabile Fumatrici espressa con valori binari. Per facilitare l’analisi se ne creerà una nuova con valori “Si” e “No”.

detach(neonati.clean)

neonati.clean = neonati.clean %>%
  mutate(Fumatrici_label = if_else(
    Fumatrici==1, "Si", "No"
  ), Fumatrici_label = as.factor(Fumatrici_label)) %>%
  dplyr::select(Anni.madre, N.gravidanze, Fumatrici, Fumatrici_label, everything())

  
attach(neonati.clean)
head(neonati.clean)
Anni.madre N.gravidanze Fumatrici Fumatrici_label Gestazione Peso Lunghezza Cranio Tipo.parto Ospedale Sesso
26 0 0 No 42 3380 490 325 Nat osp3 M
21 2 0 No 39 3150 490 345 Nat osp1 F
34 3 0 No 38 3640 500 375 Nat osp2 M
28 1 0 No 41 3690 515 365 Nat osp2 M
20 0 0 No 38 3700 480 335 Nat osp3 F
32 0 0 No 40 3200 495 340 Nat osp2 F

Verifichiamo ora che la variabile risposta Peso, oggetto dell’analisi, sia distribuita secondo una gaussiana.

Figure 8

La variabile peso sembra essere distribuita come una gaussiana, tranne che per la coda sinistra, decisamente più allungata rispetto a quella destra.

D’altronde, ci si doveva aspettare una distribuzione del genere, avendo nel campione, delle nascite avvenute con gestazioni inferiori alla norma e di conseguenza neonati con peso inferiore (cfr. Figure 2).

Conduciamo il test di Shapiro-Wilk che verifica l’aderenza dei dati ad una normale per ulteriore scrupolo.

shapiro.test(Peso)

    Shapiro-Wilk normality test

data:  Peso
W = 0.97068, p-value < 2.2e-16

Il p-value è sostanzialmente 0 e questo ci porta a rifiutare l’ipotesi nulla di normalità dei dati. Va ricordato che la variabile peso presentava comunque una forma lepticurtica e con un’asimmetria verso destra ( Table 2 ) e questo probabilmente incide sul test.

Si decide comunque di tenere a mente questo risultato ma spostare ogni valutazione una volta conseguito il modello definitivo.

4.2 Relazioni tra regressori e risposta

In questo paragrafo metteremo su grafico le correlazioni tra tutte le variabili del nostro dataset, con l’obiettivo di individuare a colpo d’occhio eventuali correlazioni che ci possano fornire una guida in sede di definizione del modello.

Figure 9: Matrice e grafico delle correlazioni per le coppie di variabili

Lasciando per ora da parte le variabili categoriche come Fumatrici, Sesso, Tipo.parto per cui questo non è il grafico migliore, per tutte le quantitative possiamo ottenere alcune informazioni:

  • Anni.madre e N.gravidanze sembrano non avere singolarmente una correlazione con il peso

  • Gestazione, lunghezza e cranio ovviamente si. In particolare sembra che la gestazione abbia un andamento non lineare rispetto al peso, con la curva che tende ad appiattirsi per alti valori. Sembra ricalcare quasi una curva logaritmica.

Per quanto riguarda invece le variabili categoriche:

Figure 10

Su tutti i boxplot possiamo notare un tratto comune: la presenza di valori molto estremi, potenziali outliars che verranno valutati in sede di modellazione. Tendenzialmente sono verso il basso e questo deve essere dovuto al fatto che nel campione è presente un certo quantitativo di parti avvenuti con gestazioni inferiori alle 37 settimane, per la precisione:

[1] 162

Quanto al sesso, abbiamo già verificato in questo paragrafo Section 3.1.1 che le differenze fisiche tra i sessi esistono e sono statisticamente rilevanti.

La variabile Fumatrici è molto sospetta perchè sembra indicare una cosa assolutamente controintuitiva, vale a dire la presenza di molti outliars sulle mamme NON fumatrici.

A un’analisi più attenta però, possiamo scartare qualsiasi considerazione fatta su questa variabile, perchè all’interno del campione abbiamo un numero di mamme fumatrici decisamente basso per poter trarre qualsiasi conclusione.

Fumatrici_label
  No   Si 
0.96 0.04 

Cerchiamo invece di capire se ci siano differenze statisticamente significative tra peso e tipo parto o ospdedale.

t.test(Peso~Tipo.parto, data = neonati.clean, alternative = "two.sided")

    Welch Two Sample t-test

data:  Peso by Tipo.parto
t = -0.13626, df = 1494.4, p-value = 0.8916
alternative hypothesis: true difference in means between group Ces and group Nat is not equal to 0
95 percent confidence interval:
 -46.44246  40.40931
sample estimates:
mean in group Ces mean in group Nat 
         3282.047          3285.063 

Il test t conferma, con un p-value a 89% che possiamo accettare l’ipotesi nulla per cui non esiste significativa differenza di peso rispetto al tipo di parto e in effetti, la cosa ha senso.

neonati.clean %>%
  filter(Ospedale != "osp3") %>%
  t.test(Peso~Ospedale, data= ., alternative = "two.sided")

    Welch Two Sample t-test

data:  Peso by Ospedale
t = -0.0092936, df = 1656.4, p-value = 0.9926
alternative hypothesis: true difference in means between group osp1 and group osp2 is not equal to 0
95 percent confidence interval:
 -51.13390  50.65161
sample estimates:
mean in group osp1 mean in group osp2 
          3270.266           3270.507 
neonati.clean %>%
  filter(Ospedale != "osp2") %>%
  t.test(Peso~Ospedale, data= ., alternative = "two.sided")

    Welch Two Sample t-test

data:  Peso by Ospedale
t = -1.6004, df = 1643.2, p-value = 0.1097
alternative hypothesis: true difference in means between group osp1 and group osp3 is not equal to 0
95 percent confidence interval:
 -92.233560   9.348156
sample estimates:
mean in group osp1 mean in group osp3 
          3270.266           3311.709 
neonati.clean %>%
  filter(Ospedale != "osp1") %>%
  t.test(Peso~Ospedale, data= ., alternative = "two.sided")

    Welch Two Sample t-test

data:  Peso by Ospedale
t = -1.623, df = 1680, p-value = 0.1048
alternative hypothesis: true difference in means between group osp2 and group osp3 is not equal to 0
95 percent confidence interval:
 -90.993927   8.590812
sample estimates:
mean in group osp2 mean in group osp3 
          3270.507           3311.709 

Anche nel confronto tra peso e ospedale il test t restituisce un p-value per tutte le 3 coppie (osp1-2, osp1-3, osp2-3) superiore al livello di significatività del 5%, ergo non rifiutiamo l’ipotesi nulla e concludiamo che la media dei pesi sugli ospedali non è diversa in maniera statisticamente significativa.

4.3 Definizione del modello con il metodo stepwise

Si partirà dal modello più complesso e si procederà per semplificazioni iterative, eliminando via via i regressori che non apportano un contenuto significativo.

4.3.1 mod 0

Partiamo dal modello che incorpora tutte le variabili, con l’unica eccezione della variabile Fumatrici che presenta la stessa informazione della variabile Fumatrici_label, con la differenza che la seconda è categorica mentre la prima è binaria.

mod0 = lm(Peso~.- Fumatrici, data = neonati.clean)
summary(mod0)

Call:
lm(formula = Peso ~ . - Fumatrici, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1123.26  -181.53   -14.45   161.05  2611.89 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -6735.7960   141.4790 -47.610  < 2e-16 ***
Anni.madre            0.8018     1.1467   0.699   0.4845    
N.gravidanze         11.3812     4.6686   2.438   0.0148 *  
Fumatrici_labelSi   -30.2741    27.5492  -1.099   0.2719    
Gestazione           32.5773     3.8208   8.526  < 2e-16 ***
Lunghezza            10.2922     0.3009  34.207  < 2e-16 ***
Cranio               10.4722     0.4263  24.567  < 2e-16 ***
Tipo.partoNat        29.6335    12.0905   2.451   0.0143 *  
Ospedaleosp2        -11.0912    13.4471  -0.825   0.4096    
Ospedaleosp3         28.2495    13.5054   2.092   0.0366 *  
SessoM               77.5723    11.1865   6.934 5.18e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274 on 2487 degrees of freedom
Multiple R-squared:  0.7289,    Adjusted R-squared:  0.7278 
F-statistic: 668.7 on 10 and 2487 DF,  p-value: < 2.2e-16

I dati salienti del modello che verranno utilizzati per tutti i confronti successivi sono:

  • il p- value dei coefficienti di regressione. In questo caso abbiamo che i coefficienti di regressioni significativamente vicini allo 0 sono:

    • gli anni della madre

    • le fumatrici

    • gli ospedali

  • R^2 adjusted pari a 0,728

Procediamo in maniera iterativa e proviamo eliminando le fumatrici, perchè effettivamente sono talmente poche rispetto al campione che non possono fornire rilievo nelle valutazioni.

4.3.2 mod 1

mod1 = update(mod0, ~.-Fumatrici_label)
summary(mod1)

Call:
lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
    Lunghezza + Cranio + Tipo.parto + Ospedale + Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1122.63  -181.96   -14.91   161.39  2615.07 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6735.4444   141.4845 -47.606  < 2e-16 ***
Anni.madre        0.8118     1.1467   0.708   0.4791    
N.gravidanze     11.1201     4.6627   2.385   0.0172 *  
Gestazione       32.3210     3.8138   8.475  < 2e-16 ***
Lunghezza        10.3064     0.3006  34.285  < 2e-16 ***
Cranio           10.4766     0.4263  24.577  < 2e-16 ***
Tipo.partoNat    29.3770    12.0888   2.430   0.0152 *  
Ospedaleosp2    -11.0363    13.4475  -0.821   0.4119    
Ospedaleosp3     28.5194    13.5038   2.112   0.0348 *  
SessoM           77.3928    11.1858   6.919 5.77e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274 on 2488 degrees of freedom
Multiple R-squared:  0.7288,    Adjusted R-squared:  0.7278 
F-statistic: 742.8 on 9 and 2488 DF,  p-value: < 2.2e-16

È rimasto tutto invariato, a conferma che effettivamente, la variabile Fumatrici non forniva un’informazione rilevante.

4.3.3 mod 2

La candidata successiva diventa la variabile Ospedale. Si sta aspettando ad eliminare quella che dai numeri sembrerebbe la prima in assoluto da scartare, gli anni della madre, per vedere se la situazione possa cambiare concentrandoci su variabili, come l’ospedale, concettualmente più separate dal peso del bambino.

mod2 = update(mod1, ~.-Ospedale)
summary(mod2)

Call:
lm(formula = Peso ~ Anni.madre + N.gravidanze + Gestazione + 
    Lunghezza + Cranio + Tipo.parto + Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1139.50  -181.60   -14.59   160.14  2633.16 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6737.9269   141.5747 -47.593  < 2e-16 ***
Anni.madre        0.8793     1.1479   0.766   0.4438    
N.gravidanze     11.4176     4.6676   2.446   0.0145 *  
Gestazione       32.6300     3.8180   8.546  < 2e-16 ***
Lunghezza        10.2839     0.3009  34.176  < 2e-16 ***
Cranio           10.4896     0.4268  24.574  < 2e-16 ***
Tipo.partoNat    30.1222    12.1038   2.489   0.0129 *  
SessoM           77.8374    11.2008   6.949 4.67e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274.4 on 2490 degrees of freedom
Multiple R-squared:  0.7278,    Adjusted R-squared:  0.727 
F-statistic: 950.9 on 7 and 2490 DF,  p-value: < 2.2e-16

Abbiamo perso circa 1 millesimo su R^2 adjusted, ora a 0,727, assolutamente accettabile e conferma definitiva che la variabile Ospedale non apportasse significativo contributo al modello.

4.3.4 mod 3

A questo punto, non resta che provare a scartare gli anni della madre e vedere che impatto otteniamo. Ce lo si aspetta contenuto, essendo quella con il coefficiente di regressione con p-value maggiore.

mod3 = update(mod2, ~.-Anni.madre)
summary(mod3)

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
    Tipo.parto + Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1129.14  -181.97   -16.26   160.95  2638.18 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   -6708.0171   136.0715 -49.298  < 2e-16 ***
N.gravidanze     12.7356     4.3385   2.935  0.00336 ** 
Gestazione       32.3253     3.7969   8.514  < 2e-16 ***
Lunghezza        10.2833     0.3009  34.177  < 2e-16 ***
Cranio           10.5063     0.4263  24.648  < 2e-16 ***
Tipo.partoNat    30.1601    12.1027   2.492  0.01277 *  
SessoM           77.9171    11.1994   6.957 4.42e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274.4 on 2491 degrees of freedom
Multiple R-squared:  0.7277,    Adjusted R-squared:  0.727 
F-statistic:  1109 on 6 and 2491 DF,  p-value: < 2.2e-16

I sospetti vengono confermati: R^2 adjusted rimane fermo a quota 0,727. Tutti i coefficienti di regressione sono ora sotto la soglia del 5% di livello di significatività, laddove il più alto p-value è a 1,2% per il tipo di parto.

4.3.5 mod 4

Proviamo ad eliminare la variabile tipo parto e vediamo cosa accade al modello. In questo caso ci aspetteremo di vedere R^2 cambiare.

mod4 = update(mod3, ~.-Tipo.parto)
summary(mod4)

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
    Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1149.37  -180.98   -15.57   163.69  2639.09 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
Cranio          10.5410     0.4265  24.717  < 2e-16 ***
SessoM          77.9807    11.2111   6.956 4.47e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274.7 on 2492 degrees of freedom
Multiple R-squared:  0.727, Adjusted R-squared:  0.7265 
F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

Si perde un ulteriore millesimo su R^2 adjusted, ora a quota 0,726. Si ritiene la perdita accettabile in cambio di una semplificazione del modello. Ad ora infatti sono state scartate 4 variabili su 8 a fronte di una perdita di 2 millesimi su R^2 adjusted.

4.3.6 mod 5

Proviamo a scartare il numero di gravidanze pregresse per valutare cosa accade al modello. È una variabile significativa, ma è quella che, pur con p-value bassissimo, lo ha ben più alto di tutte le altre.

mod5 = update(mod4, ~.-N.gravidanze)
summary(mod5)

Call:
lm(formula = Peso ~ Gestazione + Lunghezza + Cranio + Sesso, 
    data = neonati.clean)

Residuals:
    Min      1Q  Median      3Q     Max 
-1138.1  -184.4   -17.4   163.3  2626.7 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6651.6732   135.5952 -49.055  < 2e-16 ***
Gestazione     31.3262     3.7884   8.269  < 2e-16 ***
Lunghezza      10.2024     0.3009  33.909  < 2e-16 ***
Cranio         10.6706     0.4247  25.126  < 2e-16 ***
SessoM         79.1027    11.2205   7.050 2.31e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 275.1 on 2493 degrees of freedom
Multiple R-squared:  0.7261,    Adjusted R-squared:  0.7257 
F-statistic:  1652 on 4 and 2493 DF,  p-value: < 2.2e-16

R^2 adjusted rimane in effetti invariato e giustificherebbe l’eliminazione della variabile dal modello. Dal momento che gli anni della madre sono comunque ritenuti, in ambito medico, un fattore che incide sulle gravidanze, si effettueranno ulteriori test per capire se mantenere o meno la variabile.

4.3.7 Test sui modelli

Innanzitutto isoliamo le valutazioni agli ultimi 2, il modello 4 e il 5.

Il primo test che effettuiamo è l’ANOVA, acronimo di analysis of variances. Il test confronta il tasso di variazione degli scarti quadratici nel passaggio da un modello più complesso a uno più semplice, pesando questa variazione rispetto a quanto il modello più complesso riusciva a spiegare.

anova(mod5, mod4)
Res.Df RSS Df Sum of Sq F Pr(>F)
2493 188663107 NA NA NA NA
2492 188042054 1 621053.3 8.230419 0.0041541

Il p-value viene molto basso, questo porterebbe a concludere che l’ipotesi nulla per cui il modello più semplice non fa perdere informazioni rispetto a quello più complesso debba essere rifiutata. L’ANOVA dunque suggerisce di mantenere il modello 4.

Proviamo con il test di informazione bayesiano.

BIC(mod1, mod2, mod3, mod4, mod5)
df BIC
mod1 11 35208.84
mod2 9 35202.49
mod3 8 35195.25
mod4 7 35193.65
mod5 6 35194.06

Il test passa sul modello che totalizza il punteggio più basso e anche in questo caso è il modello 4 a essere preferito.

Proviamo infine a lanciare la procedura stepwise automatizzata per verificare quali passi segua e se atterri sulle nostre stesse conclusioni

step.mod = stepAIC(mod0,
        direction = "both",
        k=log(length(Peso)))
Start:  AIC=28118.61
Peso ~ (Anni.madre + N.gravidanze + Fumatrici + Fumatrici_label + 
    Gestazione + Lunghezza + Cranio + Tipo.parto + Ospedale + 
    Sesso) - Fumatrici

                  Df Sum of Sq       RSS   AIC
- Anni.madre       1     36710 186779904 28111
- Fumatrici_label  1     90677 186833870 28112
- Ospedale         2    687555 187430749 28112
- N.gravidanze     1    446244 187189438 28117
- Tipo.parto       1    451073 187194266 28117
<none>                         186743194 28119
- Sesso            1   3610705 190353899 28159
- Gestazione       1   5458852 192202046 28183
- Cranio           1  45318506 232061700 28654
- Lunghezza        1  87861708 274604902 29074

Step:  AIC=28111.28
Peso ~ N.gravidanze + Fumatrici_label + Gestazione + Lunghezza + 
    Cranio + Tipo.parto + Ospedale + Sesso

                  Df Sum of Sq       RSS   AIC
- Fumatrici_label  1     91599 186871503 28105
- Ospedale         2    693914 187473818 28105
- Tipo.parto       1    452049 187231953 28110
<none>                         186779904 28111
- N.gravidanze     1    631082 187410986 28112
+ Anni.madre       1     36710 186743194 28119
- Sesso            1   3617809 190397713 28151
- Gestazione       1   5424800 192204704 28175
- Cranio           1  45569477 232349381 28649
- Lunghezza        1  87852027 274631931 29066

Step:  AIC=28104.68
Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
    Ospedale + Sesso

                  Df Sum of Sq       RSS   AIC
- Ospedale         2    702925 187574428 28098
- Tipo.parto       1    444404 187315907 28103
<none>                         186871503 28105
- N.gravidanze     1    608136 187479640 28105
+ Fumatrici_label  1     91599 186779904 28111
+ Anni.madre       1     37633 186833870 28112
- Sesso            1   3601860 190473363 28144
- Gestazione       1   5358199 192229702 28168
- Cranio           1  45613331 232484834 28642
- Lunghezza        1  88259386 275130889 29063

Step:  AIC=28098.41
Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Tipo.parto + 
    Sesso

                  Df Sum of Sq       RSS   AIC
- Tipo.parto       1    467626 188042054 28097
<none>                         187574428 28098
- N.gravidanze     1    648873 188223301 28099
+ Ospedale         2    702925 186871503 28105
+ Fumatrici_label  1    100610 187473818 28105
+ Anni.madre       1     44184 187530244 28106
- Sesso            1   3644818 191219246 28139
- Gestazione       1   5457887 193032315 28162
- Cranio           1  45747094 233321522 28636
- Lunghezza        1  87955701 275530129 29051

Step:  AIC=28096.81
Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + Sesso

                  Df Sum of Sq       RSS   AIC
<none>                         188042054 28097
- N.gravidanze     1    621053 188663107 28097
+ Tipo.parto       1    467626 187574428 28098
+ Ospedale         2    726146 187315907 28103
+ Fumatrici_label  1     92548 187949505 28103
+ Anni.madre       1     45366 187996688 28104
- Sesso            1   3650790 191692844 28137
- Gestazione       1   5477493 193519547 28161
- Cranio           1  46098547 234140601 28637
- Lunghezza        1  87532691 275574744 29044
summary(step.mod)

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + Lunghezza + Cranio + 
    Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1149.37  -180.98   -15.57   163.69  2639.09 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -6681.7251   135.8036 -49.201  < 2e-16 ***
N.gravidanze    12.4554     4.3416   2.869  0.00415 ** 
Gestazione      32.3827     3.8008   8.520  < 2e-16 ***
Lunghezza       10.2455     0.3008  34.059  < 2e-16 ***
Cranio          10.5410     0.4265  24.717  < 2e-16 ***
SessoM          77.9807    11.2111   6.956 4.47e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274.7 on 2492 degrees of freedom
Multiple R-squared:  0.727, Adjusted R-squared:  0.7265 
F-statistic:  1327 on 5 and 2492 DF,  p-value: < 2.2e-16

Anche la procedura automatizzata atterra sul modello 4 che a questo punto teniamo come quello migliore finora.

4.3.8 Relazioni non lineari dei regressori

Osservando l’immagine Figure 9, possiamo notare come la relazione tra la variabile risposta Peso e le variabili indipendenti Gestazione, Lunghezza e Cranio, non sia proprio lineare.

Si può notare infatti un cambio di pendenza nella retta di correlazione tra le variabili, che lascia un’apertura nei confronti di relazioni delle variabili indipendenti con la variabile risposta non lineari.

Se ne indagheranno alcune per verificare se si riesca a ottenere un miglioramento del R^2 adjusted rispetto al modello 4.

4.3.8.1 Risposta logaritmica della Gestazione

L’andamento del peso rispetto alla gestazione potrebbe far intendere una relazione logaritmica. Si prova quindi ad aggiornare il modello 4 con il logaritmo della gestazione.

mod6 = lm(formula = Peso ~ N.gravidanze + I(log(Gestazione)) + Lunghezza + Cranio + 
            Sesso, data = neonati.clean)

summary(mod6)

Call:
lm(formula = Peso ~ N.gravidanze + I(log(Gestazione)) + Lunghezza + 
    Cranio + Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1153.33  -180.96   -15.52   162.10  2638.06 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -9680.4831   440.1171 -21.995  < 2e-16 ***
N.gravidanze          12.2975     4.3449   2.830  0.00469 ** 
I(log(Gestazione))  1163.2891   140.9956   8.251 2.52e-16 ***
Lunghezza             10.2554     0.3026  33.886  < 2e-16 ***
Cranio                10.5297     0.4274  24.639  < 2e-16 ***
SessoM                78.7070    11.2197   7.015 2.95e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 274.9 on 2492 degrees of freedom
Multiple R-squared:  0.7265,    Adjusted R-squared:  0.726 
F-statistic:  1324 on 5 and 2492 DF,  p-value: < 2.2e-16

Nessun miglioramento rispetto al semplice modello 4, che viene confermato anche dal BIC.

BIC(mod4, mod6)
df BIC
mod4 7 35193.65
mod6 7 35198.05

Si scarta quindi l’ipotesi di un andamento logaritmico.

4.3.8.2 Risposta quadratica rispetto alla lunghezza

Proviamo a sondare l’ipotesi che la variabile risposta sia meglio descritta con un andamento quadratico rispetto alla lunghezza.

mod7 = lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + Cranio + 
            Sesso, data = neonati.clean)

summary(mod7)

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + 
    Cranio + Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1161.44  -180.02   -11.17   165.90  2381.94 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -4.335e+03  1.442e+02 -30.059  < 2e-16 ***
N.gravidanze    1.314e+01  4.298e+00   3.057  0.00226 ** 
Gestazione      3.481e+01  3.713e+00   9.374  < 2e-16 ***
I(Lunghezza^2)  1.081e-02  3.075e-04  35.160  < 2e-16 ***
Cranio          1.047e+01  4.212e-01  24.845  < 2e-16 ***
SessoM          7.455e+01  1.110e+01   6.714 2.33e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 271.9 on 2492 degrees of freedom
Multiple R-squared:  0.7326,    Adjusted R-squared:  0.7321 
F-statistic:  1365 on 5 and 2492 DF,  p-value: < 2.2e-16

Abbiamo un incremento di ben un punto percentuale su R^2 adjusted che ora passa a 0,732, persino migliore del modello iniziale con tutte le variabili.

Verifichiamo il BIC rispetto al modello 4.

BIC(mod4, mod7)
df BIC
mod4 7 35193.65
mod7 7 35142.06

Il modello 7 scalza il modello 4 e diventa il nostro modello migliore.

4.3.8.3 Risposta quadratica rispetto al Cranio

Come fatto per la lunghezza, testiamo una relazione quadratica rispetto al Cranio, partendo però non più dal modello 4, ma dal 7.

mod8 = lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + I(Cranio^2) + 
            Sesso, data = neonati.clean)

summary(mod8)

Call:
lm(formula = Peso ~ N.gravidanze + Gestazione + I(Lunghezza^2) + 
    I(Cranio^2) + Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1157.94  -179.42   -12.32   165.84  2366.68 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)    -2.653e+03  1.176e+02 -22.570  < 2e-16 ***
N.gravidanze    1.313e+01  4.289e+00   3.062  0.00222 ** 
Gestazione      3.651e+01  3.695e+00   9.882  < 2e-16 ***
I(Lunghezza^2)  1.083e-02  3.062e-04  35.360  < 2e-16 ***
I(Cranio^2)     1.560e-02  6.218e-04  25.081  < 2e-16 ***
SessoM          7.338e+01  1.108e+01   6.620 4.38e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 271.4 on 2492 degrees of freedom
Multiple R-squared:  0.7336,    Adjusted R-squared:  0.7331 
F-statistic:  1372 on 5 and 2492 DF,  p-value: < 2.2e-16

Un miglioramento di un ulteriore millesimo di punto su R^2 adjusted che porterebbe a ritenere che il modello 8 possa sostituire il 7.

Facciamo il test BIC per averne conferma.

BIC(mod4, mod7, mod8)
df BIC
mod4 7 35193.65
mod7 7 35142.06
mod8 7 35132.63

Abbiamo un miglioramento ulteriore rispetto al modello 7 confermato anche dal BIC e pertanto teniamo il modello 8 come migliore fino ad ora.

4.3.8.4 Risposta quadratica della Gestazione

Proviamo a saggiare la possibilità di una relazione non logaritmica, ma quadratica tra la Gustazione e la variabile risposta Peso.

mod9 = lm(formula = Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) + I(Cranio^2) + Sesso, 
          data = neonati.clean)

summary(mod9)

Call:
lm(formula = Peso ~ N.gravidanze + I(Gestazione^2) + I(Lunghezza^2) + 
    I(Cranio^2) + Sesso, data = neonati.clean)

Residuals:
     Min       1Q   Median       3Q      Max 
-1156.66  -179.68   -12.29   166.35  2373.28 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.982e+03  7.042e+01 -28.145  < 2e-16 ***
N.gravidanze     1.313e+01  4.291e+00   3.059  0.00224 ** 
I(Gestazione^2)  4.840e-01  4.934e-02   9.810  < 2e-16 ***
I(Lunghezza^2)   1.087e-02  3.049e-04  35.635  < 2e-16 ***
I(Cranio^2)      1.565e-02  6.214e-04  25.179  < 2e-16 ***
SessoM           7.274e+01  1.109e+01   6.560 6.52e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 271.4 on 2492 degrees of freedom
Multiple R-squared:  0.7335,    Adjusted R-squared:  0.7329 
F-statistic:  1371 on 5 and 2492 DF,  p-value: < 2.2e-16

Vediamo che R^2 adjusted è rimasto invariato e secondo il principio del rasoio di Ockam, potremmo concludere che sia meglio mantenere il modello più semplice dal momento che non otteniamo risultati degni di nota con il quadrato della gestazione.

Effettuiamo un test bayesiano per avere ulteriore conferma

BIC(mod4, mod7, mod8, mod9)
df BIC
mod4 7 35193.65
mod7 7 35142.06
mod8 7 35132.63
mod9 7 35134.00

Il test mostra che il modello 9, che coincide con il modello 8 con l’aggiunta della risposta quadratica sulla gestazione, non è da preferirsi al modello 8.

In linea con il test bayesiano e con il principio del rasoio di Ockam, si mantiene il modello 8 come il migliore ottenuto ad ora.

4.3.9 Multicolinearità

Ora che abbiamo identificato un modello definitivo, verifichiamo che non sussista multicolinearità tra le variabili indipendenti.

Se ci fossero infatti variabili correlate tra loro, potremmo avere effetti imprevedibili sulla risposta e quindi distorsioni nelle previsioni.

vif(mod = mod8)
  N.gravidanze     Gestazione I(Lunghezza^2)    I(Cranio^2)          Sesso 
      1.023665       1.616818       1.996655       1.579492       1.041953 

La soglia comunemente utilizzata è di 5. Dal momento che tutte le variabili sono abbondantemente sotto soglia, possiamo concludere che non sussiste multicolinearità tra loro.

4.3.10 Mean Squared Error (MSE)

Verifichiamo ora che l’errore (o residuo) medio quadratico sia il minimo possibile. Il confronto lo faremo ovviamente tra il nostro modello definitivo, il modello 4 e quello iniziale.

MSE_mod0 = mean(mod0$residuals^2)
MSE_mod4 = mean(mod4$residuals^2)
MSE_mod8 = mean(mod8$residuals^2)

MSE_mod0
[1] 74757.08
MSE_mod4
[1] 75277.04
MSE_mod8
[1] 73460.61

Decisamente il modello 8 è quello da selezionare. Pur essendo più semplice, vale a dire avendo meno variabili, cosa che porta tendenzialmente a un aumento dei residui, abbiamo addirittura una riduzione rispetto al modello iniziale.

4.3.11 Commenti sul modello 8

Il modello 8 presenta dunque il giusto compromesso tra semplicità e attinenza ai dati, con relazioni sia lineari che quadratiche che nel complesso hanno portato R^2 adjusted ad essere addirittura migliore che il valore inizialmente assunto con il modello avente tutte le variabili del campione.

Questi sono i coefficienti di regressione “beta” associati ad ogni regressore

mod8$coefficients
   (Intercept)   N.gravidanze     Gestazione I(Lunghezza^2)    I(Cranio^2) 
 -2.653154e+03   1.313235e+01   3.651154e+01   1.082768e-02   1.559616e-02 
        SessoM 
  7.338031e+01 

Al netto dell’intercetta, possiamo notare come ciascuno abbia un contributo positivo sulla variabile risposta Peso. In parole povere, ogni incremento nei regressori porta un aumento nel peso.

La variabile più impattante è quella del sesso, per cui i maschi incidono con 73 grammi in più rispetto alle femmine.

Segue poi la gestazione per cui, per ogni settimana di gestazione in più, il bambino ottiene 36,5 grammi. Ovviamente questo dato è da leggersi “cum grano salis” perchè il numero di settimane di gestazione ha un suo valore naturale e non può aumentare a dismisura.

Anche il numero di gravidanze pregresse sembra avere un effetto positivo sulla risposta del Peso, laddove per ogni gravidanza in più che la mamma ha avuto, si ha un incremento di 13 grammi nel peso del bambino.

I coefficienti di regressione più piccoli sono legati alla lunghezza e al diametro del cranio. Ciò è dovuto dal fatto che queste due variabili crescono con il quadrato, e questa crescita più rapida viene bilanciata da coefficienti di regressione più bassi. In questi casi, determinare l’incremento di peso, per incremento unitario della variabile, è più complesso.

Per ogni incremento unitario di lunghezza o diametro craniale, il bambino acquisisce rispettivamente questo quantitativo di grammi in più

\[ \beta_{i}\left( 1 + 2X^{n-1}_{i} \right) \]

Dove \(X^{n-1}_{i}\) è il valore della lunghezza o diametro craniale precedente all’aumento unitario.

4.3.12 Verifica della parte erratica

Come ultimo controllo di validità del nostro modello, si deve verificare che i residui aderiscano a delle caratteristiche che li tengano “sotto controllo”:

  • distribuzione normale: i residui devono essere distribuiti normalmente, vale a dire addensati attorno allo 0

  • omoschedasticità: la varianza deve essere costante

  • non correlazione: i residui non devono essere tra loro correlati. Altrimenti si avrebbe una tendenza all’aumento o diminuzione che porterebbe il modello a divergere (come per l’omoschedasticità del resto)

Effettueremo sia una verifica grafica, che una numerica.

4.3.12.1 Correlazione dei residui

Il grafico seguente mostra, per ogni risposta del modello (asse x), quanto scarto (asse y) ci sia rispetto al corrispettivo valore reale nel dataset.

Sostanzialmente stiamo vedendo, su grafico, come siano distribuiti i residui del nostro modello. Quello che si vuole vedere è una dispersione sostanzialmente orizzontale. Non vogliamo che il modello abbia una correlazione, positiva o negativa che sia, con i residui.

Il grafico sembra confermare che non ci sia correlazione dei residui

Figure 11

Per avere una conferma numerica a quello che il grafico ci mostra, effettuiamo il test di Durbin-Watson che ha questo sistema di ipotesi. Fissiamo come sempre il livello di significatività sempre al 5%.

\[ \begin{cases} H_{0} : \mbox{Autocorrelazione dei residui}=0 \\ H_{1} : \mbox{Autocorrelazione dei residui }>0 \end{cases} \]

dwtest(mod8)

    Durbin-Watson test

data:  mod8
DW = 1.9518, p-value = 0.1139
alternative hypothesis: true autocorrelation is greater than 0

Con un p-value di poco superiore all’11% non rifiutiamo l’ipotesi nulla e confermiamo che i residui del modello non siano autocorrelati.

4.3.12.2 Distribuzione normale dei residui

Con il successivo grafico vogliamo verificare invece che i residui seguano una distribuzione normale. Per farlo utilizziamo un grafico q-q, che sta per quantile-quantile. Viene divisa in quantili la distribuzione normale e si posizionano sia sull’asse x che y. In questo modo, se i residui seguiranno questa distribuzione, dovranno posizionarsi sulla bisettrice del quadrante.

Ci si aspetta di vedere degli scostamenti, soprattutto alle code, poichè il dataset contiene bambini nati prematuramente e alcuni nati dopo le 40 settimane.

Come supposto inizialmente, i residui del modello sembrano addensarsi bene sulla bisettrice tranne che sulle code, in particolare sulla coda destra della normale (o sulla parte il alto a destra della bisettrice del grafico q-q).

Anche in questo caso effettuiamo un test numerico per vefiricare quanto i residui aderiscano a una normale standard.

shapiro.test(mod8$residuals)

    Shapiro-Wilk normality test

data:  mod8$residuals
W = 0.97856, p-value < 2.2e-16

L’ipotesi nulla del test di Shapiro-Wilk è che il vettore di dati presentato sia distribuito normalmente. Il p-value è sostanzialmente 0, quindi dobbiamo rifiutare l’ipotesi nulla: i nostri residui non aderiscono bene a una distribuzione normale standard.

Per ora teniamo a mente questo risultato e proseguiamo con l’analisi senza invalidare il modello. Tireremo a valle tutte le conclusioni.

4.3.12.3 Omoschedasticità dei residui

Terzo criterio che i nostri residui devono rispettare è quello di omoschedasticità: la varianza dei residui non deve mostre dei pattern ben precisi, ma essere sparpagliata orizzontalmente.

Sull’asse delle x troviamo le risposte del modello, mentre su quello delle y troviamo le radici quadrate dei valori assoluti dei residui “studentizzati”. Questa operazione viene fatta al fine di addensare i residui quanto più possibile per rendere più facile l’individuazione degli outliars.

Figure 12

Il grafico lascia intendere che ci sia una distribuzione dei residui abbastanza sparpagliata, ma si nota come la retta di tendenza blu sia leggermente inclinata positivamente. Effettuiamo un test numerico per eliminare ogni dubbio.

Il test di Breusch-Pagan ha questo sistema di ipotesi

\[ \begin{cases} H_{0} : \mbox{I residui sono omoshedastici } \\ H_{1} : \mbox{I residui sono eteroschedastici } \end{cases} \]

bptest(mod8)

    studentized Breusch-Pagan test

data:  mod8
BP = 58.144, df = 5, p-value = 2.938e-11

Il p-value viene bassissimo e quindi non possiamo accettare l’ipotesi nulla. Dobbiamo concludere che i residui non hanno una varianza costante, sono eteroschedastici.

Anche in questo caso, teniamo a mente il risultato e proseguiamo con l’analisi.

4.3.12.4 Leva e outliars

L’ultimo tassello da verificare è la presenza di outliars e valori di leva. Sono entrambi valori che indicano uno scostamento significativo del modello rispetto ai dati reali: nel caso degli outliars lo scostamento è sulla variabile risposta, per la leva è sui regressori.

Iniziamo con un grafico più semplice che rappresenti la distribuzione dei residui studentizzati con soglie a +2 e -2.

Figure 13

Il grafico sopra riportato mostra in effetti dei valori dei residui studentizzati oltre le soglie impostate. Questi sono punti di attenzione ma che non necessariamente si traducono in outliars perchè potrebbero non avere un effetto “dannoso” nei confronti del modello. A un’ispezione grafica infatti si riscontra come, a eccezione di pochi, siano tutti abbastanza vicini alle soglie. Verificheremo il tutto più avanti con un test specifico.

Il grafico seguente contiene molte informazioni, ma una volta chiarite dovrebbe risultare di facile lettura:

  • asse x: i valori di leva, che più sono alti e più indicano dei regressori lontani rispetto al resto

  • asse y: residui studentizzati

  • soglia leverage: linea tratteggiata blu che indica la soglia superata la quale identifichiamo i punti con leva da tenere sotto controllo. Non è detto che si traducano in outliars, ma potrebbero.

\[ \mbox{Soglia leverage}= \frac{2p}{n} \]

  • linea rossa tratteggiata: indica semplicemente lo 0 attorno al quale vogliamo si addensino i nostri residui

  • dimensione e colore dei punti: la dimensione dei punti aumenta con l’aumentare della distanza di Cook per ciascuno, come anche il colore che invece tenderà sempre più al rosso. La distanza di Cook è un parametro associato ai residui che, se superiore a 0,5 indica un punto da monitorare, se superiore a 1 indica un outliar.

Figure 14

In questo caso abbiamo valori di leva che si posizionano oltre la soglia, pochi comunque rispetto alla maggioranza, ma uno solo di loro, in alto a destra come potenziale outliar.

Effettuiamo comunque un test numerico per identificare altri outliars che possano sfuggire a un’analisi visiva del grafico

outlierTest(model = mod8)
      rstudent unadjusted p-value Bonferroni p
1549  9.028926         3.3838e-19   8.4529e-16
155   5.041006         4.9636e-07   1.2399e-03
1305  4.843831         1.3517e-06   3.3765e-03
1397 -4.287229         1.8780e-05   4.6912e-02

Vengono identificati 4 valori outliars, valori che hanno un effetto particolarmente pronunciato nei confronti del modello. Ad ogni modo, trattandosi di 4 valori su un campione di 2500 osservazioni, qualunque effetto possano avere sarà sicuramente ammortizzato dal resto del campione.

5 Previsioni del modello

Proviamo ad effettuare, alcune previsioni di test sfruttando il modello.

predict(mod8, newdata = data.frame(N.gravidanze = 2, 
                                   Gestazione = 41, 
                                   Lunghezza = 500,
                                   Cranio = 330,
                                   Sesso = "F"))
       1 
3275.426 
predict(mod8, newdata = data.frame(N.gravidanze = 0, 
                                   Gestazione = 40, 
                                   Lunghezza = 510,
                                   Cranio = 330,
                                   Sesso = "M"))
       1 
3395.389 
predict(mod8, newdata = data.frame(N.gravidanze = 1, 
                                   Gestazione = 40, 
                                   Lunghezza = 500,
                                   Cranio = 330,
                                   Sesso = "M"))
       1 
3299.162 
predict(mod8, newdata = data.frame(N.gravidanze = 1, 
                                   Gestazione = 40, 
                                   Lunghezza = 480,
                                   Cranio = 330,
                                   Sesso = "M"))
      1 
3086.94 

I test sono stati condotti su nascite di cui si conosce il peso e, tranne nel primo caso, in cui la bambina è nata pesando 2,8 kg, quasi mezzo chilo di differenza rispetto al modello, negli altri casi la corrispondenza è molto vicina alla realtà.

Proviamo infine una previsione con dati parziali. Supponiamo ad esempio di voler prevedere il peso di una bambina la cui madre è alla terza gravidanza e partorirà alla 39ma settimana.

In questo caso, mancando i dati di lunghezza e cranio, utilizzeremo le medie del nostro campione.

predict(mod8, newdata = data.frame(N.gravidanze = 2, 
                                   Gestazione = 39, 
                                   Lunghezza = mean(neonati.clean$Lunghezza),
                                   Cranio = mean(neonati.clean$Cranio),
                                   Sesso = "F"))
       1 
3250.079 

6 Conclusioni

Possiamo concludere che il modello si adatta bene sui valori “normali” e nascite non particolari, mentre sulle code, parti prematuri o tardivi, o anche pesi eccessivamente alti o bassi pur in gestazioni normali, non garantisce ottime performance.

Le evidenze riscontrate sull’eteroschedasticità e sulla distribuzione normale dei residui si consiglia di non farle incidere sulla valutazione del modello. Effettivamente, è il migliore identificabile, rimanendo nella linearità, con i dati a disposizione.

Si suggerisce pertanto di avviare una fase di “valutazione sul campo” del modello, usandolo nei vari reparti dei 3 ospedali per fare previsioni sul peso alla nascita. Il peso effettivo verrà poi confrontato con la previsione e verranno tratte nuove conclusioni sulle performance a valle.

Sarà inoltre possibile, in questo modo, fare una valutazione dei benefici apportati in termini economici tramite l’efficientamento delle risorse. Se i benefici apportati saranno comunque superiori dei costi sostenuti per i casi di previsione errata, sarà comunque accettabile mantenere un modello di cui si conoscono esattamente i limiti.